﻿#light

module DwarfFSharpMpi

open System
open System.Collections.Generic
open MPI
open System.Runtime.InteropServices
open DwarfConfigurator

[<DllImport(@"CliTools.dll")>]
extern void start()

[<DllImport(@"CliTools.dll")>]
extern double finish()

let processGraphModels(model : DataModel) : DataModel =
    //temporary storage for max value for psi array
    let mutable  maxvalind = 0   
    //temporary storage for max value for psi array                     
    let mutable maxval = 0.0
    let mutable tempVal = 0.0 
    // 1. Initialization  delta and psi 
    let numberOfStates = model.hmm.getNumberOfStates()
    
    //only for root thread
    if (model.isRootThread) then
        for i = 0 to  numberOfStates - 1 do
            model.vit.delta.[0].[i] <- model.hmm.initialStates.[i] * (model.hmm.observationProbabilityMatrix.[i].[model.vit.observationSequence.[0]])
            model.vit.psi.[0].[i] <- 0
            
    // 2. Compute 
    //sending to all instances the initial delta and psi values
    Communicator.world.Broadcast<int>(ref model.vit.psi.[0], DataModel.ROOT_PROCESS)
    Communicator.world.Broadcast<double>(ref model.vit.delta.[0], DataModel.ROOT_PROCESS)

   
    let lengthOfObsSeq = model.vit.getLengthOfObservationSequence()
    
    let workDeltaArray = Array.create model.lengthForProc 0.0  //temporary storage of delta values
    let workPsiArray = Array.create model.lengthForProc 0      //temporary storage of psi values
    let mutable p = 0                                          //helper index
    
    for t = 1 to lengthOfObsSeq - 1 do
        //wait for all processors
        p <- 0
        for j = model.startIndx to model.endIndx - 1 do
            maxval <- 0.0
            maxvalind <- 0
            //finding the maximum
            for i = 0 to numberOfStates - 1 do
                tempVal <- (model.vit.delta.[t - 1].[i] * (model.hmm.stateTransitionMatrix.[i].[j]))
                if (tempVal > maxval) then
                    maxval <- tempVal
                    maxvalind <- i
            workDeltaArray.[p] <- maxval * (model.hmm.observationProbabilityMatrix.[j].[model.vit.observationSequence.[t]])
            workPsiArray.[p] <- maxvalind
            p <- p + 1
        //send to all instances the results of computation on step t
        Communicator.world.GatherFlattened(workPsiArray,model.arrayOflengthForProc,DataModel.ROOT_PROCESS,ref model.vit.psi.[t])
        Communicator.world.GatherFlattened(workDeltaArray, model.arrayOflengthForProc, DataModel.ROOT_PROCESS, ref model.vit.delta.[t])
        Communicator.world.Broadcast<double>(ref model.vit.delta.[t], DataModel.ROOT_PROCESS)

    // 3. Termination
    //final computation on master process
    if (model.isRootThread) then
        model.vit.setProbability(0.0);

        model.vit.sequenceOfState.[lengthOfObsSeq - 1] <- 0
        for  i = 0 to numberOfStates - 1 do
            if (model.vit.delta.[lengthOfObsSeq - 1].[i] > model.vit.getProbability()) then
                model.vit.setProbability(model.vit.delta.[lengthOfObsSeq - 1].[i])
                model.vit.sequenceOfState.[lengthOfObsSeq - 1] <- i
  
        // 4. Path (state sequence) backtracking
        let backtrackingArray = [lengthOfObsSeq - 2 .. -1 .. 0]
        for t in backtrackingArray do
            model.vit.sequenceOfState.[t] <- model.vit.psi.[t + 1].[model.vit.sequenceOfState.[t + 1]]

    model
  
let main() =
    let args = System.Environment.GetCommandLineArgs()
    using(new Environment(ref args))(fun _->
        let conf = new Configurator()
        let mutable model = conf.configure(args,true)
        if(model.isRootThread) then
            model.DWARF_NAME <- "Graph Models, fsharp Mpi kernel."
            model.DWARF_NAME_OUTPUT <- "Dwarf.FSharp.Mpi"
            start()
            
        model <- processGraphModels(model)
        if(model.isRootThread) then
            let time = finish()
            conf.Finish(model,time)
    )
  
main()   